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Very accurate wave functions are calculated for small transition metal oxide molecules. These 
wave functions are decomposed using reduced density matrices to study the underlying correlation 
of electrons. The correlation is primarily of left-right type between the transition metals and the 
oxygen atoms, which is mediated by excitations from the nominal single Slater ground state into 
antibonding and d-type orbitals. In a localized representation, this correlation manifests itself in a 
2-electron hopping term that is off-diagonal. This term is of similar magnitude to the commonly 
considered Hubbard-type on-site interaction. 



One of the grand challenges in modern condensed 
matter theory is the description and prediction of the 
properties of correlated electrons. Of particular inter- 
est are transition metal oxides, which exhibit effects in- 
cluding high Tc superconductivity p4 |2j, Mott insulator 
behavior[3], and colossal magnetoresistance[ll IHl, all of 
which owe their existence to electron correlation effects. 
Control of these strong correlation effects has the po- 
tential to open up many new areas in both physics re- 
search and technology, much like the control of weakly 
correlated electronic structure has enabled innumerable 
advances in the past 80 years since the development of 
that theory. Much of the research in physics to date 
has concentrated on the development of phenomenolog- 
ical models of strong correlation, such as the Hubbard 
model, that have had many successes in helping to un- 
derstand these systems. However, when considered from 
first principles, the underlying Coulomb Hamiltonian of 
a strongly correlated system is the same as for a weakly 
correlated system. The difference between the two is an 
emergent property of the many-electron wave function. 

There has been a large amount of effort devoted to 
treating strongly correlated systems starting from the 
first principles Hamiltonian. These efforts have ranged 
from phenomenological corrections to density functional 
theory (DFT), such as DFT+U[i [7]/DFT-hDMFT[H], 
to other extensions of DFT using hybrid functionals, to 
GW perturbation theory 9., quantum chemistry [lOj, and 
quantum Monte CarlopTMlS). While all these approaches 
have had varying levels of success, there still remains a 
gap: very few calculations have been performed that ap- 
proach the exact solution to the Schrodinger equation on 
strongly correlated systems and analyze the nature of the 
correlated wave function in these challenging materials. 

This article is meant as a first step at understand- 
ing the relevant correlations in transition metal oxides. 
Quantum Monte Carlo(QMC) is chosen as a vehicle to 
do this because of two major considerations. First, the 
strong dynamic correlation that is present in transition 
metal oxides is easily described using explicit correla- 



tion, which is efficiently evaluated using Monte Carlo 
techniques. Second, and perhaps more importantly, the 
quantum Monte Carlo methods are able to perform cal- 
culations on extended systems efficiently, which is unique 
for an explicitly correlated wave function based method. 
Learning what elements of the wave function are nec- 
essary for accurate treatment of transition metal oxide 
molecules in QMC with thus provide a valuable guide for 
larger molecules and extended systems. In this article, 
we will explore in what ways electron correlation breaks 
single particle symmetry and thus discover what terms 
should be in an effective model of electron correlation. 



I. SYMMETRY AND ELECTRON 
CORRELATION 

Let's start witha discussion of how one particle symme- 
tries are broken with electron correlation. Suppose the 
one-particle Hamiltonian Hn, has a symmetry such that 
it commutes with the single particle symmetry operators 
Ai and _Bj, where i refers to the single-particle electron 
number. We will consider without loss of generality only 
one symmetry operator A. Then the eigenstates of Hn, 
can be labeled based on their symmetries as follows 

\s) = \E;ai,a2,a3,a4,...), (1) 

where E is the eigenvalue of Hib, and {a} is the set of 
one-electron eigenvalues of A. 

Now suppose that we add to the Hamiltonian a two- 
particle effective interaction H2b so that the full Hamil- 
tonian is Hii, + H2b- Further suppose that H2b commutes 
with the single-particle operator Ai. Then for two eigen- 
states of Hib \si) and \sj), we can label them as 

|s,;) = |i?„{a},), (2) 

Then, 

{■3^\H2b\sJ) ^0, (3) 

if {a}, ^ {a}, (4) 
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That is, i?2b does not change the block-diagonal nature 
of Hib if it commutes with Ai. This can occur when H2b 



2 



is diagonal in any basis that preserves the symmetry of 
Hib, among other cases. 

Let's now contrast a reahstic real-space interacting 
Hamiltonian with the above H2b- For electrons, the ab- 
initio Hamiltonian is Hn, + '^^j ^/fij in atomic units. 
The l/vij term conserves the overall symmetry of the 
system, but does not conserve particle-by-particle sym- 
metries. More explicitly, the eigenstates can be labeled 
with the total symmetry of the state X^i^ii where the 
sum is over all electrons, but not the individual Ai's. 

For example, if the system has cylindrical symmetry, 
the total angular momentum is a good quantum number 
for the interacting Hamiltonian, but the angular momen- 
tum of a particular electron is not. If one attempts to em- 
ulate the effect of the l/rij term by using an interaction 
term that conserves the single particle angular momen- 
tum, then one is enforcing the symmetry on a particle- 
by-particle basis. Because conserving the one-particle 
symmetry aids in solving the model system, many com- 
monly used models for electron correlation obey Eqn [3] 
for at least some single particle symmetries. 

For modeling transition metal oxides, a very com- 
mon effective interaction is the on-site d-orbital H = 
'^iUnJnj, where hi is the number operator on an 
atomic-like d-orbital. This interaction explicitly does 
not allow eigenstates that are mixtures of single parti- 
cle rotational states. For example, if the nominal ground 
state of the transition metal monoxide MnO is (3d^2p'^, 
2p'^), where the states before/after the commas indi- 
cate spin up/down, then superimposing the configura- 
tion (3d^4s^,2p'^,2p^3d^), which involves a double elec- 
tron hopping, is not allowed. We shall see from accurate 
calculations of the electronic structure of first-principles 
systems that this superposition is critical when consider- 
ing the first principles Hamiltonian. 



II. METHOD 

To obtain accurate first-principles results, we use 
variational quantum Monte Carlo (VMC) and fixed- 
node diffusion Monte Carlo^ (FN-DMC). Quantum 
Monte Carlo methods are well-described elsewhere in the 
literature jl5], so they will be described very briefly here. 
VMC is a straightforward implementation of the varia- 
tional method using Monte Carlo evaluation of the en- 
ergy expectation value. FN-DMC simulates the imagi- 
nary time Schrodinger equation is simulated to obtain the 
lowest energy state consistent with a given nodal surface, 
which further improves over the variational results. 

The core electrons are replaced with a 
pseudopotential[Tni HZ]) which is treated in the lo- 
cality approximation [T8] . In the transition metals, the 3s 
and 3p electrons are considered part of the valence. All 
calculations are performed using the QWalk[T^ package. 
The trial wave function is the multi Slater Jastrow(MSJ) 



wave function: 



^{K) = exp(C/) J2 CfeDet[0f >(r,)], 



(5) 



where the determinants to include are taken from a 
configuration interaction in singles and doubles calcula- 
tion, the one-particle orbitals are taken from a hybrid 
B3LYP 20, calculation in GAMESS[ni[22], and the Jas- 
trow factor U is the one described in Ref [TT]. The 
coefhcicnts Ck are energy optimized !23j simultaneously 
with the Jastrow parameters. A similar approach has 
been shown to efficiently produce high accuracy on a 
benchmark set of molecules. Enough determinants were 
included that the one-particle density matrices did not 
change upon including more determinants. 



A. Calculation of the reduced density matrices 

The reduced density matrices are evaluated in quan- 
tum Monte Carlo using the following integrals for the 
single particle reduced density matrix (1-RDM): 

f^l = E / <l>l{r'a)MraW{K)^iR)dr'JR (6) 

and for the two particle reduced density matrix (2-RDM) : 

[ ^Ur'aWM)Ura)Mra) (7) 

where R = (ri,r2, . . . ,rAr), R'^ = (ri,r2, . . . ,r^, . . . ,rAr), 
and i?"^ = (ri, r2, . . . , r^, . . . , r^, . . . , rjv) and normaliza- 
tion is omitted. These matrices can be spin resolved, re- 
sulting in two 1-RDMs for spin up and down, and three 
2-RDMs for the combination of up/up, up/down, and 
down/ down. 

In VMC, these can be evaluated by sampling two ad- 
ditional coordinates and in addition to the many- 
electron coordinate R. In this work, the additional 
coordinates are sampled from the distribution /(r) = 
^■0f(r). R is drawn as usual from The ex- 

pectation values are thus given as follows (after a some 
rearrangement of terms and inclusion of normalization) : 



N,M. 



(8) 



and similarly: 



Pij,kl ~ Ar Ar Ar a7 ' ^) 



with 



N,, = 




(10) 
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TPSSh[25] 


3.48 


3.43 


3.58 


3.97 




RMC(SJ) 


4.61(5) 


4.11(5) 


4.64(5) 


4.76(4) 


5.3(1) 


DMC(MSJ) 


3.77(2) 


3.16(2) 


3.89(5) 


3.27(4) 


4.92(4) 


Exp[2a 


4.55 


3.34(1) ;28; 


3.355 


3.88 





TABLE I: Dipole moments in Debye. The fixed-node RMC re- 
sults have been obtained with a single deteriminant of B3LYP 
orbitals. 



Choosing /(r) properly to sample and rj, increases 
the efhciency substantially, as well as using symmetry to 
evaluate ^'(i?a) for all a for a Slater- Jastrow or multi 
Slater-Jastrow wave function with little work. The full 
implementation can be found in the QWalk code. The 
density matrices are expressed in a basis of B3LYP one- 
particle orbitals, except where noted in the text. The 
density matrices using the mixed estimator in DMC are 
indistinguishable from the VMC results when for the con- 
verged wave functions in this work, which further rein- 
forces the accuracy of the wave functions. 

The above procedure was performed for the early tran- 
sition metal monoxides ScO,TiO,VO, CrO, and MnO, 
and the late transition metal dioxides Mn02, Fe02, and 
C0O2. The latter set has an interesting transition from a 
bent to straight bond that is very sensitive to the treat- 
ment of electron correlation. 



III. RESULTS AND DISCUSSION 

A. Geometry and dipole moments 

For the transition metal monoxides, the dipole mo- 
ment (Table |l]) is very challenging to calculate using quan- 
tum Monte Carlo. For the accurate wave functions con- 
sidered here, the agreement with experiment is much bet- 
ter than using a Slater-Jastrow wave function, although 
it still appears quite difficult to converge the dipole mo- 
ment, since the energy of a state is not very sensitive to 
the dipole moment. 

In the case of the transition metal oxides, the bond an- 
gle (Table [Tl| is predicted poorly [29 by most GGA meth- 
ods and even hybrid methods, so this is a stringent test 
of the treatment of electron correlation. Diffusion Monte 
Carlo with a Slater-Jastrow nodal surface underestimates 
the bond angle of Fe02 , and gives a flat potential energy 
surface for C0O2 around the linear geometry. A more ac- 
curate wave function fixes these defects and clearly agrees 
with the bond angles obtained in experiment. 
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TABLE II: Minimum energy bond angles at r=1.6Afor the 
transition metal dioxide molecules. The DMC methods both 
obtain the correct bond length of 1.6 A. 



B. Occupation numbers of the one-particle density 
matrix 

In the basis of B3LYP orbitals, the single particle re- 
duced density matrix is very accurately diagonal; there is 
likely little to gain in this case in orbital optimization of a 
single Slater determinant, beyond using B3LYP orbitals. 
Thus we only report the diagonal elements (Fig [T]) . In 
the case of the monoxides, the a-symmetry orbitals have 
the lowest occupation number in the nominally occupied 
set of states, while in the virtual space, the up electrons 
occupy mostly the a* orbital and the down electrons oc- 
cupy a number of virtual orbitals. Interestingly, the d-like 
singly occupied orbitals have occupation numbers closer 
to 1 than the bonding-type orbitals for the spin majority, 
indicating that in some sense, these orbitals once occu- 
pied are not that strongly correlated. 

CrO is a special case because it has a degenerate 
ground state in the single particle approximation. Cor- 
relation lifts this degeneracy by mixing the two states, 
which is responsible for some of the outliers in Fig[l] 

C. Off-diagonal elements of the two particle 
reduced density matrix 

To analyze the breaking of one-particle symmetry by 
the interaction, we can turn to the two particle reduced 
density matrix. Suppose that we expand the a state \'^) 
in terms of a basis of Slater determinants. Since this 
matrix can be written in second quantized form in a basis 
as 

pgL = (*i44c.c,i*)> (11) 

one can show that if |5')'s expansion contains two Slater 
determinants \s) and |s') such that 

W)^clclc.cj\s), (12) 

then pfjf.^ is nonzero if (i,j) 7^ {k,£). Conversely, if 
there are no Slater determinants in the expansion of 
connected by Eqn |12[ then the matrix element is zero. 
We will use this to detect the satisfaction or lack thereof 
in Eqn [3] for different symmetry classes. 
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FIG. 1: (color online) The 1-RDM for transition metal monox- 
ides in the B3LYP Kohn-Sham basis. Stochastic errors are 
smaller than the symbol sizes. The diagonals only are plot- 
ted, since the off-diagonal elements are very small and do not 
change the picture. 



Most of the non-zero off-diagonal elements of the two- 
particle reduced density matrix involve the cr-bonding- 
like orbitals, labeled in Fig [2] The only exception here 
is again CrO with the degenerate single-particle ground 
state. For all of the systems considered here, the spin- 
like off-diagonal elements are larger than the spin-unlike, 
which is surprising-perturbation theory implies that the 
spin unlike correlation should be larger. For most of the 
materials, the off-diagonal elements are in the form of an 
exchange between the cr-like orbitals and another sym- 
metry orbital, which breaks the one-particle rotational 
symmetry. These elements are of the form of a 2-electron 
hopping; a typical example of which would be 



(13) 
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where cti and (T2 are two different cr-symmetry states and 
TT is a TT-symmetry state. 

One sees a clear trend in the monoxide molecules; 
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FIG. 2: The first column is tti the second is and the 
third is fi- All molecules are at their equilibrium geometries, 
except C0O2, which is set to 140° bond angle for easier sym- 
metry comparison. 



as the state goes from a doublet (ScO) to a sextuplet 
(MnO), the correlation increases monotonically. In the 
dioxides, as the state goes from a quadruplet (Mn02) to a 
doublet (C0O2), the correlation decreases monotonically. 
This is likely the reason for the increase in bond angle 
through this series, since there are fewer empty states in 
C0O2 in which to perform exchanges, the electrons re- 
pulse each other more as the bond angle closes, therefore 
tending towards a 180 degree bond. This can be seen in 
the occupation number of one of the main virtual orbitals 
(Fig [3]) as a function of angle. 

The picture emerging from the calculations can be 
summarized as follows. The cr-like bond between the 
transition metal and the oxygen experiences a strong dy- 
namic interaction with other electrons. This bond is thus 
most likely to be partially occupied. In terms of virtual 
excitations, the most likely excitation is an electron ex- 
cited from a singly occupied d-like state into a low-lying 
virtual orbital, and then an electron occupied from the cr- 
like orbital into the newly de-occupied d-like state. Inter- 
estingly, this scenario cannot be described with an on-site 
Hubbard U-like term, due to the symmetry, as discussed 
above. We can see the effects of this when trying to fit a 
low-energy model to the physics of the monoxide MnO. 



D. Fitting a model: MnO 

It is interesting to consider the minimal effective model 
that can reproduce the density matrices considered in 
this paper. As an example, the case of MnO is consid- 
ered, with the valence space made from the Mn 3d and 
4s states and the oxygen 2p states. Including the 4s state 
is essential to reproduce the physics, since the partial oc- 
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FIG. 3: (top) The occupation number of the 17th orbital 
(summed over spin) for Mn02, Fe02, and C0O2 as a function 
of angle at a bond length of I.60A. The lines are guides to 
the eye. (bottom) The 17th orbital for Mn02, with a positive 
isosurface in blue, and negative in red. The other TM-O2 
molecules are qualitatively similar. 



cupation of this state is large. The full state space in this 
case is only 756 states, so it can be solved by exact diag- 
onalization, and the model parameters can be optimized 
to reproduce the 2-RDM diagonals. The Hamiltonian 
considered is 



H = E + T + Ui, 



2-cxchangc ? 



(14) 



where E is the one-particle energy of the localized 
orbitals (3 parameters for 2p,3d,and 4s), T is the 
one-particle hopping parameter (one parameter for tt- 
symmetry orbitals, and 3 parameters for cr-symmetry or- 
bitals) , ?7intrasito IS the Hubbard U and onsitc interactions 
that are diagonal in the local basis, and J/2-cxchango is a 2- 
electron hopping term that is off-diagonal in the localized 
basis: 



2-oxchango — t/Cl C2 '^2 ''Cs ^Ci i 

C1C2 



(15) 



where C is a spin and site index {i,a). More details on the 
precise Hamiltonian is available in the Appendix. The 
last term is critical; without it, the 2-RDM is very poorly 
reproduced, having RMS errors of ± 0.17 on numbers 
that vary between zero and one. With the last term, the 



RMS errors are reduced to ± 0.01. It thus appears that a 
2-electron exchange term is critical to describe the corre- 
lation between the transition metal and the oxygen. This 
intersite exchange term appears to be very rarely con- 
sidered in theoretical descriptions of strongly correlated 
materials, usually entering only in an intrasitc form [30], 
although it has been noted [31' that a similar term can 
result from downfolding a Hubbard model with intersite 
Coulomb interaction to a t-J model. 



IV. CONCLUSION 

Given accurate many-body wave functions for small 
transition metal oxide molecules, the largest correlations 
break single-particle rotational symmetry. In an effec- 
tive model that reproduces the two-body physics of the 
MnO molecule, this effect is similar to the size of the 
on-site Hubbard-like interaction. It appears that to ac- 
curately describe the electron interactions in these ma- 
terials, while a Hubbard-like U term can aid in obtain- 
ing rough agreement with the true ground state, accu- 
rate agreement requires breaking the one-body rotational 
symmetry. It remains to be seen whether or not these ef- 
fects are more or less important as the system size grows 
larger. This is under current investigation. 

If it is true that the correlations presented here are 
generally important in transition metal oxide systems, 
then they may provide a guide to building more accu- 
rate trial wave functions to use in quantum Monte Carlo 
calculations. The basic two-particle hopping could be de- 
scribed with a compact wave function of only a relatively 
few Slater determinants, and the optimized with power- 
ful techniques. Investigations on this front are also under 
development. 

The methods presented here are quite general, and can 
be applied to both solids and larger molecules. The 1- 
RDM is very inexpensive to evaluate with a proper imple- 
mentation, while the 2-RDM is somewhat more expen- 
sive, but can be made to scale as the number of electrons 
squared if localized basis functions are used. The key 
concept here is to use the reduced density matrices to 
carve the large Hilbert space into pieces that are more 
easily analyzed, and to combine this with the accurate 
wave functions attainable using quantum Monte Carlo 
techniques. Rather than relying only on energetics to fit 
models, the information provided by a single calculation 
can inform models of the electron correlation to increase 
the physical realism. 
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V. APPENDIX 



A. Effective multiband model for the MnO 
molecule 



Orbital indices: 



7 






2p cr 


1 


Mn 3d (7 


2 


Mn 4s cr 


3 


2p TT^ 


4 


Mn 3d -Kx 


5 


2p TTy 


6 


Mn 3d TTy 


7 


Mn 3d ^1 


8 


Mn 3d 52 



Fitted parameters for MnO 



E2p 


-14.97 




-15.63 


Eis 


-8.85 




-4.87 


Toi 


-4.41 


To2 


-5.1 


Tl2 


0.75 


U2p 


-0.32 


Um 


4.02 




5.11 


U2p-2p 


0.27 


Usd-Sd 


0.09 


Usd-is 


-1.32 


C^lOOl 


-2.65 


^^2002 


-2.19 


C^2112 


4.96 



H — E + T + J/intrasite + U2 -exchange ) 

(16) 



E= E2pn,+ i;3dn«+ XI ^4sn» (17) 

ie{0,3,5} ie{l,4,6,7,8} ie{2} 



f= Y T^clcj+Toiclci (18) 

(i,j)e{(3,4),(5,6)} 

+ To2clc2 + Ti2c\c2 + h.c. 



t^hubbard= ^2pninl+ Y Uadnjn} (19) 

ie{0,3,5} iG{l,4,6,7,8} 

+ Y 

ie{2} 

(20) 

C^intrasite = Yj C^2p-2p^i' flj (21) 

i>je{0,3,5},(7,cr' 

+ Y Usd-3dnfn^ (22) 

»,je{l,4,6,7,8},<T,<7' 

+ Y Uzd-isfil hj + h.c. 

j6{2},je{l,4,6,7,8},CT,<T' 

(23) 



U2 -exchange — 

Uloolc[clcoCi + U2002clclcoC2 + U2112clc[ciC2 

(24) 



